What is claimed is: 



1 . A method of simulating multibody dynamics of a molecular system so as to produce 
output data representative of a time evolution of the molecular system, comprising: 

providing a set of equations for characterizing multibody dynamics of the molecular 
system, the set of equations being constrained by a constraining equation; 

providing an integrator for integrating the set of equations, wherein the integrator is 
tailored for the set of equations so as to satisfy the constraining equation; and, 

applying the integrator to the set of equations over a time step, so as to produce output 
data representative of the time evolution of the molecular system relative to the time step. 

2. A method according to claim 1, further including applying the integrator to the set of 
equations over a predetermined number of time steps, so as to produce output data representative 
of a time evolution of the molecular system relative to a time interval corresponding to the 
predetermined number of time steps. 

3. A method according to claim 1, wherein the constraining equation includes |e| = 1 . 
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4. A method according to claim 1 , wherein the set of equations for characterizing multibody 
dynamics includes equations of motion for each of a plurality of bodies in the system, the 
equations of motion given by 

e = W(u)) e 
x = C(e)u 

Pu, = Gjq) + p T [A„]p 
Pu = G u (q) + p T [A u ]p 
P; = Gt(q) + p T [Az]p 

wherein the matrices W(co) and C(e) represent the kinematical relations and are written 
in the form 



IV (co) = 1 
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and, 



C(e) = 



1-2(ei + eS) 
2(e 2 ei-e 3 e 0 ) 
2(e 3 ei + e2e 0 ) 



2(eie 2 + e 3 e 0 ) 
1-2(ei + e?) 
2(e 3 e 2 -eie 0 ) 



2(eie 3 -e 2 e 0 ) 
2(e 2 e 3 + eie 0 ) 
1-2( ei 2 + e 2 ) 



, respectively. 



5 . A method according to claim 1 , wherein applying the integrator to the set of equations 
over a time step further includes propagating Euler parameters of the set of equations for 
characterizing multibody dynamics, according to the matrix exponential of e m = E (cq V2 ) e 0 ■ 
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6. A method according to claim 1, wherein applying the integrator to the set of equations 
over a time step further includes: 

propagating a momentum vector to a half step f = Af /2 , while freezing one or more 
position states at f = 0 ; 

propagating a vector of rigid position variables to the half step f = Af /2, while freezing 
one or more flexible position variables at f = 0 , and the momentum vector at f = Af /2 ; 

propagating the vector of flexible position variables to the full step t = At using momenta, 
recomputed in modal velocities, at f = At/2 : 

propagating the vector of rigid position variables to the full step t = Af , while freezing the 
flexible position variables at t = Af and the momentum vector at f = Af /2 ; and, 

propagating the momentum vector to the full step f = Af freezing the position states at 

f = Af. 

7. A method according to claim 1, wherein applying the integrator to the set of equations 
over a time step further includes decomposing the set of equations for characterizing multibody 
dynamics first with respect to momentum states, and then with respect to position and rigid and 
flexible degrees of freedom. 

8. A method according to claim 1, wherein applying the integrator to the set of equations 
over a time step further includes decomposing the set of equations for characterizing multibody 
dynamics first with respect to rigid and flexible degrees of freedom, and then with respect to 
position and momentum states. 
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9. A method according to claim 1, wherein providing an integrator includes solving a 
kinematic equation associated with the set of equations for characterizing multibody dynamics 
via a solution based on a matrix exponential of the form 

E (a> V2 ) = e wM f = cos (yAt) I + sin (yAt ) D . 

10. A computer system for simulating multibody dynamics of a molecular system so as to 
produce output data representative of a time evolution of the molecular system, comprising: 

an analytical structure for characterizing multibody dynamics of the molecular system, 
wherein the analytical structure is constrained by a constraining equation; 

an integrator for integrating the analytical structure, wherein the integrator is tailored for 
the structure so as to satisfy the constraining equation; 

wherein the integrator integrates the structure over a time step, so as to produce output 
data representative of a time evolution of the molecular system relative to the time step. 

11. A computer system according to claim 1 0, wherein the analytical structure includes a set 
of equations modeled in a computer code. 
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12. A computer system according to claim 11, wherein the set of equations includes 
equations of motion for each of a plurality of bodies in the molecular system, the equations of 
motion given by 

e = W(uj) e 
x = C(e)u 

P„ = Gjq) + p T [Aa,]p 
Pu = Gu(q) + p T [Au]p 
Pz = Gt(Q) + P T Mp 

wherein the matrices W(co) and C(e) represent the kinematical relations and are written 
in the form 



l/V(u>)4 
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and, 



C(e) = 



1 - 2 (el + ei) 2 ( ei e 2 + e 3 e 0 ) 2 (e^a - e 2 e 0 ) 
2(e 2 ei - e 3 e 0 ) 1 - 2 (ei + e?) 2 ( e2 e 3 + eie 0 ) 
2(e 3 ei + e 2 e 0 ) 2(e 3 e 2 -eie 0 ) 1 -2( e ? + ei) 



respectively. 



13. A computer system according to claim 10, wherein the constraining equation includes 
Iel = 1. 



14. A computer system according to claim 10, wherein the integrator further integrates the 
structure over a predetermined number of time step, so as to produce output data representative 
of a time evolution of the molecular system relative to a time interval corresponding to the 
predetermined number of time steps. 
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15. A computer system according to claim 10, wherein the integrator further propagates Euler 
parameters of the set of equations for characterizing multibody dynamics, according to the matrix 
exponential of evz = E (o)v 2 ) e 0 • 

16. A computer system according to claim 10, wherein the integrator further solves a 
kinematic equation associated with the set of equations for characterizing multibody dynamics 
via a solution based on a matrix exponential of the form 

E (^i /2 ) = e w(6,V2) f = cos (yAt ) / + sin (yAt) D . 

17. A computer system according to claim 10, wherein the integrator further decomposes the 
set of equations for characterizing multibody dynamics first with respect to momentum states, 
and then with respect to position and rigid and flexible degrees of freedom. 

18. A computer system according to claim 10, wherein the integrator further decomposes the 
set of equations for characterizing multibody dynamics first with respect to rigid and flexible 
degrees of freedom, and then with respect to position and momentum states. 
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19. A computer system for simulating multibody dynamics of a molecular system so as to 
produce output data representative of a time evolution of the molecular system, comprising: 

means for characterizing multibody dynamics of the molecular system, wherein the 
analytical structure is constrained by a constraining equation; 

means for integrating the analytical structure, wherein the integrator is tailored for the 
structure so as to satisfy the constraining equation; 

wherein the means for integrating integrates the structure over a time step, so as to 
produce output data representative of a time evolution of the molecular system relative to the 
time step. 
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20. A computer system according to claim 19, wherein the analytical structure includes 
equations of motion, modeled in computer code, for each of a plurality of bodies in the molecular 
system, the equations of motion given by 



t=u t 

P„ = Gjq) + p T [A a )p 
Pu = Gu(q) + p T [Au]p 
P 4 = Gt(q) + p T [A 4 ~]p 

wherein the matrices W(co) and C(e) represent the kinematical relations and are written 
in the form 



e = W(io) e 
x = C(e)u 
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and, 



Wy " OJz 



Wz - co y 



1-2(el + e3) 2( ei e2 + e 3 eo) 2( ei e3-e 2 eo) 
C (e) = 2 (e 2 ei - e 3 e 0 ) 1 - 2 ( e | + e?) 2 (e 2 e 3 + eie 0 ) , respectively. 
2(e3ei + e 2 eo) 2(e 3 e 2 -eie 0 ) 1 -2( e ? + el) 
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